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ABSTRACT 

We construct three dimensional axisymmetric, cuspy density distributions, whose po- 
tentials are of Stackel form in parabolic coordinates. As in Sridhar and Touma (1997), 
a black hole of arbitrary mass may be added at the centre, without destroying the 
Stackel form of the potentials. The construction uses a classic method, originally due 
to Kuzmin (1956), which is here extended to parabolic coordinates. The models are 
highly oblate, and the cusps are "weak", with the density, p oc l/r fc , where < k < 1 . 

Key words: galaxies: elliptical and lenticular, cD — galaxies: kinematic and 
dynamics — galaxies: structure 
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1 INTRODUCTION 

Kuzmin(1956) showed that the density at a general point 
of an oblate, axisymmetric mass model, whose gravitational 
potential is of Stackel form in elliptic coordinates, is sim- 
ply related to the density on the axis of symmetry. Later, 
Kuzmin (1973) generalized the theorem to the fully tri- 
axial case. A corollary of Kuzmin's theorem is that the 
density so constructed is automatically non-negative ev- 
erywhere as long as it is non-negative on the short axis. 
The interest in Stackel models arises from the fact that 
the motion of point masses in such potentials is completely 
integrable by the Hamilton-Jacobi method; these poten- 
tials are also termed "separable", since the dynamics in 
three dimensions separates into three independent oscilla- 
tions (c.f. de Zeeuw 1985a, b for a clear discussion of orbits, 
and Kuzmin's theorem). Constructing distribution func- 
tions for models of galaxies in a steady state (c.f. Bin- 
ney & Tremaine 1987) is facilitated by the existence of 
three isolating integrals for three dimensional models. Re- 
cently, Sridhar and Touma (1997) — hereafter referred to as 
ST — presented a family of cuspy, scale-free mass models 
of non-axisymmetric discs whose potentials are of Stackel 
form in parabolic coordinates. The surface density behaves 
like E oc 1/r 7 , where < 7 < 1. Here we generalize 
Kuzmin (1956) to parabolic coordinates and naturally ex- 
tend the non-axisymmetric discs of ST to three dimensional, 
axisymmetric, cuspy density distributions, whose potentials 
are of Stackel form in parabolic coordinates. Our main result 
follows in § 2. Scale-free cusps are presented in § 3, and a 
complete classification of orbit families is given in § 4. A dis- 
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cussion of the results, 
models. 



in § 5, completes this paper on mass 



2 KUZMIN'S CONSTRUCTION APPLIED TO 
PARABOLIC COORDINATES 

For axisymmetric configurations, the mass density p, and 
the gravitational potenti al U, th at it generates are both 
functions only of R = -\/x 2 + y 2 and z, the cylindrical ra- 
dius and the coordinate along the axis of symmetry, respec- 
tively. Since the potential is independent of the azimuthal 
angle ((/>), the component of the angular momentum par- 
allel to the z-axis, L z , is conserved. Hence the dynamics 
may be described as two-dimensional motion in a merid- 
ional (<f> = constant) plane, under the influence of an effec- 
tive potential defined by 



I/off = U(R, z) + 



2R 2 ' 



(1) 



We study cases when the motion in the meridional plane 
is separable in parabolic coordinates, (f, rf). These may be 
defined as the roots for r of the equation 



t 2 - 2zt - R 2 = 0. 



(2) 



Solving this quadratic equation gives 
£ = z + y/z 2 + R 2 

77 = z~y / z 2 + R 2 , (3) 
where we have chosen r\ < < £. The inverse relations are 



R 2 = -£77, 



(4) 
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where r = \/R 2 + z 2 is the spherical radius. Surfaces of con- 
stant £ and 77 are paraboloids of revolution which cut the 
2-axis at £/2 and 77/2, and the x-y plane at R = £ and R — 77 
respectively. These surfaces are orthogonal to the meridional 
plane and mutually orthogonal as well; thus (£, 77, 4>) form an 
orthogonal coordinate system. Positive and negative ranges 
of the 2-axis are covered by £ and r\ respectively: 



£ = 0, 2 = 2r], for 2 < 
T) = 0, 2 = 2£, for 2 > 0. 



(5) 



The most general axisymmctric potential in parabolic 
coordinates, for which the Hamilton-Jacobi equation sepa- 
rates, has the Stackel form, given by 



U(£,v) = 



f+(0 , Fjn) 
£-v v~£' 



(6) 



If U(£,r)) is of Stackel form, so is the effective potential, 

U c fi, which is defined in equation (|l|).^] The potential due to 
a point mass is obtained by choosing F- — F+ = 2GM. It is 
important to note that, to describe an interesting case like 
the Kepler problem, F+ 7^ F_ ; in general, we let F+ and F_ 
be quite different functions of their arguments. The density 
corresponding to the point mass assumes the form, 



M5(r) = M 



\ «(£-v) J 



Let p and U be the density and potential, respectively, of a 
smooth distribution of matter^]. If we now include a point 
mass at the origin, 



Ptot 
Utot 



p + p. 



U 



2GM 

£-v' 



(8) 



are the total density and potential. These are related to each 
other through Poisson's equation, or equivalently 



V 2 U = 4nGp, 



(9) 



where the contribution of the point mass has been dropped 
from both sides. Expressed in parabolic coordinates, Pois- 
son's equation for (the axisymmctric) potential and density 
reads 



£-v 



d_ d_ 



&_ d_ 



(10) 



We substitute the expression given in equation (JSJ) for U, 
in equation ( |lo| ) . This results in the following expression for 
the smoothly distributed density of matter: 



nGp(£,V) = 



£F+ + r)F_ 

{£ - n? 



£ + v 
(£-vY 



(F + ~F_). 



(11) 



where the primes denote differentiation with respect to the 
arguments, £ or 77, as the case may be. It proves convenient 
to employ the single variable r in place of the two variables 



£ and rj. Then the density is completely determined by 
F+(£) ifr>0 



Ft) 



F-(rj) ifr<0 



(12) 



which is evidently one function of one variable. It is also 
useful to define constants C T and D T , that depend only on 
the behaviour of F for small Irj: 



C T 



D T = 



C+ = lim 7) _ r]F_ if r > 
C- = lim 5 _ £F+ if t < 

D+ = \im v ^ (vF-)' if r > 
D- =]im e _ >0 (|F+)' ifr<0 



(13) 



(14) 



We are now ready to apply Kuzmin (1956) 's method. 
Let us assume that the density profile on the 2-axis, ip(r), 
is specified: 



i>{ T ) = 



i>+(O = P (£,0) ifr>0 
•tp-(n) = p(0,T)) ifr<0 



(15) 



Then we can use equation ( pd| ) to first determine F(t), and 
hence p({, 77). Restricting equation (|ll|) to the 2-axis, results 
in the following (ordinary) differential equation for F(t) : 



F - — = nGrip(T) - ^ - — , (((») 

T T T 

where primes denote differentiation with respect to r. We 
define ^(r) = 4>(s) ds, and integrate equation ( fl6| ) twice, 
to obtain the general solution, 

F(t)=tvG / s<S>(s)ds + A T + B T — + ^ln|r|+D T r,(17) 
Jo 2 2 

where A T and B T axe constants of integration. We recall, 
from the discussion following equation (|6j), that A v — A% 
will merely add to M; to describe a smooth mass distribu- 
tion, we may, without loss of generality, set A„ — A^ — . 
Substituting this solution for F(t) in equation (hlh provides 
us with a general expression for the Stackel density (whose 
2-axis profile is ip(r)) in parabolic coordinates: 



p(£,v) 



e^+(£)+v 2 4>-(v) 



(£ - r,f 



2£V 



+ 



(£-v) 3 

Pbcd ■ 



-*-(!?)] 



(18) 

We recall that the Stackel nature of the potential — hence 
separability of the Hamilton-Jacobi equation — remains un- 
affected if a point mass (see eqns. jjj] and with arbitrary 
M is added at the origin. In equation ([18^ Pbcd is the con- 
tribution from the B T , C T , and D T terms in equation (|l7|). 
We show in the Appendix that this contribution is generally 
unphysical, so we ignore pbcd in the rest of the paper. It 
is evident that the remaining expression for p(£, 77) is non 
negative. 



t To verify this we only need note that adding L 2 /2^ and L 2 /2r] 
to F+ and F— respectively, adds the centrifugal term to U(^,tj) . 
§ We include density cusps in our "smooth" distribution! 



3 SCALE-FREE CUSPS 

A 2-axis density profile of the form ip(r) = po (?"o/M) fe , can 
give rise to a finite, non-negative Stackel density only if 
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< k < 1; this is so because we require cuspy densities to 
have k > 0, and only when k < 1 is 



*(t) 



(19) 



finite. Here, 6(t) is ±1 accordingly as r is positive, or neg- 
ative. F(t) is readily computed using equation (0). Scale- 
free cusps, with a central black hole are described by 



-2K\r)\ : 



GM 
k + GM 



(20) 

where K = ttGpqTq /2(1 — k)(2 — k). The potential-density 
pairs have simple expressions in parabolic coordinates, 



2K[( i - k + \ v \ :i - k ] 



2GM 



Ps(£,v) = Por 



fcC +\m 



(t-v) 2 



PoW) 2£r) r i_fc 
(1 - fc) (f - T7) 3 L« + ™ J 



(1-fc) (C-r?) 3 
+A75(r) , (21) 

or equally well in spherical polar coordinates: 

U s {r,9) = Kr 2 - k [{l + cosej s - k + {l-cos6) 3 - k ]- — 



p B {r,6) = Po[^) [(2-fc-cos6>)(l 

+ (2-fc + cos#)(l-cos6>) 5 
+MS(v) . 



(22) 



Meridional sections of the isocontours of the volume den- 
sity and potential for k equal to 0.1 and 0.5 are shown in 
Figure 1. For definiteness, the section may be taken to be 
the x-z plane. The density (Figures la and lc) is stacked 
on highly oblate figures of revolution, that display a dip 
near the z axis. In all figures, the black hole is shown as a 
dot at the centre. In Figures lb and Id, the contribution 
of the black hole to the potentials is suppressed. As in the 
corresponding figures of ST, the potential isocontours have 
fc-independent axis ratio of 2. 

The motion of point masses in all axisymmetric po- 
tential respects two classical, isolating integrals of motion: 
(E,L Z ), the energy and the angular momentum about the 
z-axis. Stackel potentials are special in that the Hamilton- 
Jacobi equation separates in special, orthogonal coordinates, 
giving rise to a third isolating integral I3 (c.f. Landau & Lif- 
shitz 1976). For motion in the meridional plane, L z may be 
treated as a given constant. We give below, without deriva- 
tion, expressions for E and I3: 



E = 



h 



2£ 



2 1 

Pi + 



_2rj_ 



(23) 



2r 1 p 2 l + F-(r ] ) + ^--Er ! . 



(24) 







Figure 1. Isocontours of Density and Potential: Figures (a) and 
(c) show isocontours of the surface density for two different values 
of 7. The corresponding potentials are displayed in the panels on 
the right in (b) and (d). Successive isocontours of the density 
have ratios of 1.1 and 1.5 in (a) and (c) respectively. Successive 
potential isocontours in (b) and (d) have ratios of 2. The location 
of the central black hole is shown as a solid dot in all four figures, 
although the contribution to the potentials is not included in (b) 
and (d). 



4 ORBITS 

Stackel potentials being integrable, orbits are confined to 3- 
tori in the 6-dimensional phase space. The three isolating 
integrals, E, L z and 73, allow us to classify orbit families. 
When the black hole is present, E can take all real values, 
whereas only positive values are permited for cusps without 
black holesP| For given E, the magnitude of L z must be less 
than, or equal to the angular momentum of the circular, 
equatorial orbit of the same energy. Having fixed E and L z , 
the third integral, 1%, will determine the excursions in £ and 



77. From equation (u 
non negative gives 



u), the requirement that and p v be 



h > 9(0, 



-h>g(\v\)- 



where 



2Ks 3 - k + L 



2s 



Es - GM : 



s > 0. 



(25) 



(26) 



Thus, the range of 73 is determined by the minimum of the 
function g(s), which is necessarily non positive. If this min- 
imum value is denoted by — 7 m (75, L 2 ), we obtain the con- 
dition — 7 m < 73 < 7 m . We note that the dynamics is truly 
scale-free only in the absence of the black hole. 

Given (E, L z ,I-i), an orbit belonging to this 3-torus can 
explore a region of real space, whose boundaries are deter- 

^ No fundamental significance is attached to these positive values 
of E; they merely reflect the fact that our cuspy potentials have 
been chosen to be non negative. 



© 0000 RAS, MNRAS 000, [|-| 



mined by the intersections of coordinate surfaces, x-z cuts 
of these regions, for representative orbit families, are shown 
shaded in Figures 2a-2d. For definiteness, we have chosen a 
scale-free potential of the form given in equation (^) . When 
^3 flips sign, the allowed region for the orbit is reflected 
about the x-z plane; hence we will consider only I3 > . 

1. The L z 7^ orbits are all loops about the 2-axis. These 
avoid the 2-axis, and generally fill a region of space that is 
rotationally symmetric about the z-axis; Figure 2a displays 
the intersection of such a region, for a generic loop, with the 
x-z plane. When I3 = 0, the region is obviously symmet- 
ric about the x-y plane, as is evident from Figure 2b. The 
shaded region now is filled with resonant orbits; the two 
frequencies, for motion in the meridional plane, are equal. 

2. Polar orbits, obtained when L z = 0, display greater vari- 
ety. Motion is restricted to a plane which, for convenience, 
we think of as the x-z plane, and the orbits can cross the 
z-axis. In this plane, the orbits are similar to those that live 
in the phase space of the non-axisymmetric disc potentials 
considered by ST. To keep the presentation self-contained, 
we go over the classification here: 

(i) E < : g(s) increases monotonically with s, reaching a 
minimum value of —GM for s = 0; thus —GM < I3 < GM. 
For a fixed (say, positive) value of 73, the motion is bounded 
by the coordinate curves 

< £ < Cm , 

< \v\ < Vm (27) 

where £ m (E,Iz) > rj m {E,I^) > are the two roots of 
g(s) — ±13. The orbits, named lenses by ST, fill a lenticular 
region bounded by the parabolas £ = £ m and r\ — —r) m - As 
the shaded region in Figure 2c indicates, the lens orbits visit 
the origin. When | 1 = GM, its maximum value, the shaded 
region collapses to an interval of the z-axis, < z < £ m /2. 
For negative values of J3, the roles of £ and r\ are inter- 
changed. 

(ii) E > 0: g(s) is not monotonic in s; a minimum, — J m < 0, 
is attained at s = [E/2K(3 - k)] 1/{2 ~ k) . Lenses occur for 
|/3 1 < GM, but a new type of orbits appears when GM < 
I ^3 1 < Im' the orbits are centrophobic bananas (they cross 
the z-axis, but avoid the origin), and the allowed region 
is shown in Figure 2d. They are parented by a resonant 
banana orbit which is obtained when I3 = J m . When M. = 
0, E is always positive; the lens orbits disappear, leaving the 
bananas as the only generic family of orbits. 



5 DISCUSSION 

We have generalized Kuzmin's theorem to parabolic coor- 
dinates and constructed three dimensional, axisymmetric, 
cuspy mass models whose potentials are of Stackel form. A 
point mass may be added at the centre, while leaving the 
separable (i.e. Stackel) nature of the potentials unimpaired. 
We have also presented a set of scale-free cusps, p oc l/r k , 
where < k < 1, corresponds to the weak cusps seen in the 
"core-type" elliptical galaxies of Lauer et. al (1995). 

Kuzmin's original construction in elliptic coordinates 
allows the central axis ratio to be specified as an indepen- 
dent parameter, whereas the very choice of parabolic coor- 
dinates forces the isopotential contours (when the contri- 
bution of the black hole to the potential is neglected) to 



Figure 2 




Figure 2. Orbit Families: Intersection of the region of real space 
filled by orbits on a 3-torus (in phase space) with the x-z plane is 
shown as the shaded area in all the figures, (a) and (b) correspond 
to loop orbits which have L z ^ 0. (a) is for a "generic" set of loops, 
each of which densely fills the shaded region, (b) corresponds to 
loops with J3 = ; these resonant, simple, closed (1:1) orbits. 
Both (c) and (d) are polar orbits with L z = 0. These live entirely 
on one meridional plane, here taken to be the x-z plane, (c) and 
(d) are filled by centrophilic lenses, and centrophobic bananas, 
quite similar to those discussed in ST. 



be significantly oblate, with axis ratio equal to 2. There are 
some other essential differences between models in parabolic 
coordinates on the one hand, and those studied earlier 
in elliptic/ellipsoidal coordinates. As de Zeeuw, Peletier & 
Franx (1986) noted, Stackel models in ellipsoidal coordinates 
do not give rise to cuspy densities; furthermore, the Stackel 
nature is completely destroyed if a point mass is included at 
the centre. Our models are fortunate to escape these limita- 
tions. However, the axisymmetric models we have been able 
to construct have densities that can be no steeper than 1/r; 
this not only limits our models to weak cusps, but does not 
allow construction of finite mass models. On the other hand, 
the Stackel models in ellipsoidal coordinates have densities 
that fall off as steeply as 1/r 4 far from the centre. 

When the potential of the black hole is neglected, the 
potential isocontours of our models are similar, concentric 
ovals. They may be compared with the potentials of the 
power-law galaxies of Evans (1994), which are stratified on 
similar, concentric spheroids. However, the dynamics in the 
power-law models is not integrable (there is no third in- 
tegral). And like all scale-free models they have infinite 
mass. Evans' oblate, power-law galaxies also have isodensity 
contours that are dimpled near the z-axis (see also Rich- 
stone 1980, and Toomre 1982 for discussions of scale-free 
models with flat rotation curves; with k — 2 for these mod- 
els, the cusps are very strong indeed). Evans notes that there 
is no real evidence that the density distribution of ellipticals 
is not dimpled, and further adds that his strongly dimpled 
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models have projected surface density isocontours that are 
somewhat boxy. Since our models have only weak cusps, the 
surface density can be computed only with the specifica- 
tion of a truncation at large r. While it may be tempting 
to speculate that our Stackel models capture some essential 
features of the dynamics in the nuclear regions of giant el- 
lipticals (whose isophotes are known to be boxy), we must 
be cautious: truncation at large r will likely spoil the exact 
separability that our models presently enjoy. Whether one 
can effect a truncation that retains much of the integrable 
nature is yet to be explored. 

While three dimensional, axisymmetric cusps are not 
expected to be chaotic (c.f. Richstone 1982), it is never- 
theless useful to have explicit integrable mass models with 
simple analytic expressions (for the density and potential) 
that also include a central black hole. The loops are the 
only generic orbit family, a fact that was recognised by Ol- 
longren (1962) in his classic study of orbits in galactic po- 
tentials, [jj] Polar orbits fall in two categories: the centropho- 
bic bananas and the centrophilic lenses. The models pre- 
sented in this paper are the natural extension of the non- 
axisymmetric two dimensional cusps that we have recently 
presented in ST. Together they exhaust the family of scale- 
free cusps that can be constructed in parabolic coordinates. 
Integrable, triaxial cusps still pose a challenge. 
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APPENDIX A: 

We describe the contribution to p of the B T , CV and D T 
terms in Eq. (|l7|). Writing pbed = Pb + pc + pd, the three 
contributions in order are: 

B+_B_ ( -2fr? 



Pb(g,t?)= " (t^Ts )■ ( A1 ) 



which, in the usual (r, d, (f>) polar coordinates, takes the form 
, m f B +-B-\ sin 2 6> 

The density pb is positive for B + > B- and zero on the 
z-axis, as is expected of a homogeneous solution of Eq. (|l6|). 

KM = (^)^- I ^[^Ml + co S , ) 
-C_ln(l-cos0)], (A3) 
, n\ ( D- — D + \ cos8 . 

The density pa changes sign across the x-y plane, unless we 
require that D+ or D- be discontinuous. All terms in pbed 
are unphysical and deserve to be ignored. 



II Having determined the boundary of the region (of the merid- 
ional plane) in which the orbit is confined, Ollongren (1962) loses 
interest in parabolic coordinates as being useful for a description 
of integrable stellar motion in our Galaxy. 
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